Astronomy & Astrophysics manuscript no. KruipPCI2010 


©ESO 2010 


February 25, 2010 





Mathematical properties of the Simplex algorithm 

C. J. H. Kruip, 1 J.-P. Paardekooper 1 B. J. F. Clauwens, & V. Icke 1 '* 

'Sterrewacht Leiden, Postbus 9513, 2300RA Leiden, The Netherlands 
e-mail: kruip@strw.leidenuniv.nl 

Received *** 2009; accepted *** 

ABSTRACT 

Context. Analytical and numerical analysis of the Simplex radiative transfer algorithm, which features transport on a Delaunay 
triangulation. 

Aims. Verify whether the SimpleX radiative transfer algorithm conforms to mathematical expectations, to develop error analysis and 
present improvements upon earlier versions of the code. 
Methods. Voronoi-Delaunay tessellation, classical Markov theory. 

Results. Quantitative description of the error properties of the SimpleX method. Numerical validation of the method and verification 
of the analytical results. Improvements in accuracy and speed of the method. 

Conclusions. It is possible to transport particles such as photons in a physically correct manner with the SimpleX algorithm. This 
requires the use of weighting schemes or the modification of the point process underlying the transport grid. We have explored and 
applied several possibilities. 
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1. Introduction 

At present one of the major challenges in computational astro- 
physics is to correctly account for radiative transfer in realistic 
macroscopic simulations. Whether one considers the formation 
of single stars or the evolution of merging galaxies, incorpora- 
tion of radiative transfer is the next essential step of physical re- 
alism needed for a deeper understanding of the underlying mech- 
anisms. 

With the advent of the computer as an important catalyst, the 
last century has seen a myriad of efforts to solve the equations of 
radiative transport within a numerical framework. The resulting 
algorithms cover a wide range of applications where generally a 
method is tailored to a specific problem. In most cases, speciali- 
sation either means the choice for a specific physical scale (con- 
sider detailed models of stellar atmospheres or large-scale cos- 
mological simulations) or emphasis on physical processes rele- 
vant for the problem at hand. 

1.1. The SimpleX algorithm 

The SimpleX algorithm for radiative transfer, which is the sub- 
ject of this text, has been designed to steer clear from these limi- 
tations. Due to a modular structure, different physical processes 
can be included straightforwardly allowing different areas of ap- 
plication. More importantly, the method has no inherent refer- 
ence to physical scale and thus can be applied to problems with 
typical length scales that lie many orders of magnitude apart. 
Moreover, due to its adaptive nature, many orders of magni- 
tude in optical depth and spatial resolution can be resolved in 
the same simulation. 

Conceived by Ritzerveld & Icke (2006) and implemented 
by Ritzerveld ( 2007} , the SimpleX algorithm solves the general 
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equations of particle transport by expressing them as a walk on 
a graph. In this sense the method can be thought of as a Markov 
Chain on a closed graph. Transported quantities travel from node 
to node on the graph, where each transition has a given proba- 
bility. 

This approach has several advantages which we would like 
to highlight here. More specifically SimpleX 

- does not increase its computational effort or memory use 
with the number of sources in a simulation and consequently 
treats for instance scattering by dust and diffuse recombina- 
tion radiation without added computational effort 

- naturally adapts its resolution to capture the relevant physical 
scales (e.g. expressed in photon mean free path lengths) 

- works in parallel on distributed memory machines 

- is compatible with grid-based as well as point based hydro- 
dynamics codes (where the latter is the more natural com- 
bination due to the point-based nature of both SimpleX and 
SPH) 

- is computationally cheap due to the local nature of the 
Delaunay grid 

1.2. Error analysis 

Covering many orders of magnitude in spatial resolution (and 
optical depth) is a considerable challenge for radiative transfer 
methods in general. Often, approximations need to be employed 
to keep the problem tractable, in turn making a robust error anal- 
ysis difficult. 

For the vast majority of numerical methods the errors are 
measured by comparison with a fiducial run of the code. This 
is usually a simulation wherein many more time steps and/or a 
higher spatial resolution is used than would be normally feasible. 
Usually, some sort of convergence to a solution is observed and 
this solution is accepted as the correct one, at least within the 
possibilities of the method. 
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For SimpleX we cannot do such a convergence test in gen- 
eraQ When the spatial resolution is increased far beyond the 
resolution dictated by the local mean free path length of the pho- 
tons, several effects to be described in Sec. [3] will tend to make 
the radiation field more diffuse, decreasing instead of increasing 
accuracy. This property of SimpleX is not a weakness but instead 
inherent to the fact that the method uses a physical grid wherein 
deviations from its natural resolution often imply deterioration 
of the solution. 

A great advantage is that the SimpleX algorithm, due to its 
mathematically transparent nature, allows for an analytical as- 
sessment of its error properties. The resulting prescriptions are 
quite general and can be applied to different regions in parame- 
ter space, an advantage over the 'numerical converge approach' 
usually applied in the error analysis of radiative transfer meth- 
ods. The development of analytical descriptions and discussion 
of their implications will be the main focus of this text. 

1.3. Outline 

Given the desirable properties of the SimpleX algorithm, we aim 
to assess if the various inaccuracies that arise inevitably in a nu- 
merical method can counterweigh SimpleXs virtues in any case. 

We will show that the use of high dynamic range Delaunay 
triangulations as the basis for radiative transfer introduces sys- 
tematic errors which reveal themselves as four distinct effects: 
diffusive drift, diffusive clustering, ballistic deflection and bal- 
listic decollimation. 

We will describe and quantify these effects and, subse- 
quently, provide suitable solutions and strategies. 

Following this, we will demonstrate how the derived mea- 
sures can be used to constrain the translation of a physical prob- 
lem to a transport graph in such a way as to avoid or minimise 
errors. 

Finally, we will discuss the results in a broader picture and 
outline related and future work; specifically the recent imple- 
mentation of a dynamically updating grid in SimpleX. 

Although developed in the context of SimpleX, many of our 
results are relevant for other transport algorithms that work on 
non-uniform grids (e.g. AMR grids or SPH particle sets.) 



2. Radiative transfer on unstructured grids 

In this section we describe the three means of transport that con- 
stitute the SimpleX algorithm. Starting from a general descrip- 
tion of the equation of radiative transfer in inhomogeneous me- 
dia, we will delve into the specifics of SimpleX. This in turn asks 
for an introduction to the Delaunay triangulation which lies at 
the heart of our method and plays a central role in this text. 



2.1. General radiation transfer 

The equation of radiative transfer for a medium whose properties 
can change both in space and time is given by 
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1 Although this is possible if the direction conserving transport to be 
introduced in Sec. |2. 5. 3] is used. Alternatively, several instances of the 
same simulation with a different random seed for the grid construction 
can be averaged to obtain error estimates as well. 



Here I v = I(t, x, Ji, v) is the monochromatic specific intensity, h 
is a unit vector along the ray and rj v and Xv are the monochro- 
matic emission and opacity coefficients respectively. 

If - ^ <k 1 (in other words: if I v is not explicitly time de- 
pendent or the time and space discretisation is such that c can be 
considered infinite) this equation simplifies to 
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If we take the spatial derivative along the ray and divide by Xv, 
Eq.Q can be rewritten as 
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where we have defined the source function S v = t] v /x v and op- 
tical depth dr = Xvds and s parametrises the distance along the 
ray. Eq.Q can be solved numerically if rj v andxv are known lo- 
cally. What we mean by 'local' depends on the type of discreti- 
sation of the volume. The form of Eq.([3]l suggests that the optical 
depth is the most natural variable to discretise. Doing this will 
in general result in cells of unequal volume and, hence, an irreg- 
ular computational mesh. As mentioned in Sec. [T] and explained 
in more detail in Sec.[6]we generate a set of nuclei which reflects 
the underying optical depth field and connect these nuclei with a 
Delaunay triangulation. 

2.2. A natural scale 

The fundamental idea behind the SimpleX algorithm is the fact 
that there exist a natural scale for the description of radiative 
processes: the photons local mean free path 
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Here <j is the total cross section and n the number density of the 
particles responsible for the extinction. As photons travel typ- 
ically one mean free path before interacting with the medium, 
information on a much smaller scale does not necessarily yield 
a deeper understanding of the physical problem at hand. 

Following this line of reasoning, the next step is to choose a 
computational mesh that inherently carries this natural scale. In 
other words: use an irregular grid which resolution adapts locally 
to the mean free path of the photons traveling over 



2.3. The grid 
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This computational mesh or transport graph, is constructed by 
first defining a point process that represents the underlying 
(physical) problem and, second, by connecting these points ac- 
cording to a suitable prescription. 

Specifically, the point process is defined by prescribing the 
local point field density as a function of the scattering and ab- 
sorption properties of the matter through which the particles 
propagate. The resulting points are connected by means of the 
unique Delaunay triangulation ( |Delone|1934 i. Thus, the points 
that carry the triangulation, which we call nuclei in the remain- 
der of this text, are connected in a graph, along whose connecting 
lines (called edges) the particles are required to travel. 

This choice for constructing a transport graph has several 
advantages: first, the connection with the physical processes is 



2 In our current implementation we use a single frequency bin and 
an effective cross section so the mean free path depends on the number 
density and the point-to-point distance only. 
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evident; second, the Delaunay triangulation (and its dual, the 
Voronoi tessellation ( Voronoi 1908 1) is unique and is in many re- 



spects the optimal unstructured partitioning of space (Schaap & 
|Van de We ygaert 2000); third, the Voronoi-Delaunay structures 
carry all relevant information needed for the transport process, 
making the transport step itself very efficient. 

We note that the use of Voronoi-Delaunay grids does not im- 
ply that the transfer of photons has to proceed along the edges of 
the Delaunay triangulation. Another approach would be to trace 
rays through either the Voronoi cells or Delaunay simplices. This 
ray-tracing could be done classically with long-characteristics 
(e.g | Mihalas & Weibel Mihalas|| 19841 ) or using Monte Carlo 

As said, in SimpleX, the 



Abbott & Lucy][19851 l. 



methods (e.g 

photons travel from cell to cell interacting with the medium as 
they go along. In this sense, the transport of radiation is treated 
as a local phenomenon and the global nature of the radiative 
transfer problem is dealt with by sufficient iterations of this local 
transport. The intimate relation between the structure of the grid 
and the transport of the photons is one of the reasons why the 
SimpleX algorithm is exceptionally fast. 

The construction of the grid itself is a task performed by ded- 
icated software which is freely available on the web. As said, 
once given the generating nuclei, the Voronoi-Delaunay struc- 
ture is unique. It is in creating the distribution of generating nu- 
clei that we can manipulate the properties of our computational 
grid. The translation from a given density (opacity) field either 
as particles or as cells to a Voronoi-Delaynay grid optimal for 
SimpleX is a fundamental part of the algorithm. 
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Fig. 1. Voronoi tessellation of the plane (solid lines). Each cell 
contains all points that are closer to its nucleus (indicated by 
a dot) than to any other nucleus. The corresponding Delaunay 
triangulation is shown in dashed lines. Note: only visible nuclei 
are included in the triangulation. 



2.4. Voronoi-Delaunay structures 

Because it plays a central role in this text, we now proceed with 
a concise introduction to the Voronoi-Delaunay grid. 

Let a set of nuclei be given in a D-dimensional space with a 
distance measure (in all of our applications, we use the isotropic 
'Pythagoras' measure). Partition this space by assigning every 
one of its points to its nearest nucleus. All the points in space 
assigned to a particular nucleus n form the Voronoi region of 
nucleus n (if D - 2, the region is often called Voronoi tile). 
By definition, the Voronoi region of n consists of all points of 
the given space that are at least as close to n (according to the 
distance measure of that space) as to any of the other nuclei. 
The set of all points that have more than one nearest nucleus 
constitute the Voronoi diagram of the set of nuclei (Edelsbrunner 
2006 ; |Okabe et al.||2000l |O'Rourke||200Tl |Vedel Jensen|1998) . 
For clarity, most of our explanatory illustrations will use D = 2 
(see Fig.[T]for an example in the plane), but our applications have 
D = 3. 

The set of all points that have exactly two nearest nuclei 
n\ , ri2 is the Voronoi wall between these nuclei. If D — 2, such 
a wall is a line; with D = 3 it is a plane; and so on. In our 
astrophysical applications, D = 3, in which case (pathological 
configurations excepted) triples of planar walls come together 
in lines, the points of which have three nearest nuclei; the lines, 
finally, join in quadruples at nodes, single points that have four 
nearest nuclei. A Voronoi wall thus separates the Voronoi regions 
around two nuclei. The connection between these nuclei, called 
an edge, is the geometric dual of the wall. The set of all such 
edges is the Delaunay triangulation ( |Delone|1934[ ) of the point 
set. The Delaunay triangulation can be defined as the unique tri- 
angulation of nuclei for which the interior of a Delaunay simplex 
(triangle in two dimensions) contains no other nuclei. 

In our algorithm, no other connections between nuclei are al- 
lowed, hence only adjacent Voronoi regions are connected. The 



expectation value for the number of these Delaunay neighbours, 
A, is 6 in 2D and 15.54- • ■ in 3D. Voronoi regions based on an 
isotropic distance measure are convex; for a point process con- 
taining nuclei, the number of geometric entities (walls, edges, 
nodes, etc.) is O(N). 

Most important for our application to transport theory, is that 
the Delaunay triangulation has a minimax property: this trian- 
gulation is the one that has the largest smallest angle between 
adjacent triangle edges. That is to say, of all possible triangula- 
tions of a given point set, the Delaunay triangulation is the one 
that maximises the expectation value of the smallest angle of its 
triangles. Colloquially: the Delaunay triangulation has the least 
'sliver-like' triangles, and the most 'fat' triangles. 

Another advantageous property of the Delaunay triangula- 
tion is that every edge pierces the wall between the Voronoi re- 
gions it connects at right angles. This justifies the use of Eq.Q 
where the spatial derivative is taken along the ray. 

2.5. Three types of transport 

In SimpleX, transport is always between neighbouring Voronoi 
cells (they are the ones connected by the Delaunay triangula- 
tion). On average every nucleus has 6 neighbours in 2D and 
15.54- • • neighbours in 3D. When photons travel through a cell, 
the optical path length, /, is taken to be the average length of the 
Delaunay edges of that cell. If the number density of atoms in 
the cell is given by n the fraction of photons that are removed 
from the bundle, A^ rem , is given by 

JV rem = N in e-'^ (5) 

where <x is the total cross section which may be due to multiple 
extinction processes. 

In general, extinction can be subdivided into absorption 
and scattering. Radiation that is removed by absorption pro- 
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cesses will in general change the temperature of the medium and 
change its chemical state but it need not be transported further at 
this point. 

Photons that are removed from the bundle due to scattering 
should propagate to neighbouring cells, either isotropically or 
with a certain directionality^] This can be accomplished using 
the diffuse transport method schematically depicted in the left 
panel of Fig.[2jand explained in more detail in Sec. |2. 5.1 



The radiation that is not removed from the bundle, however, 
needs to travel straight on along the original incoming direction. 
For optically thick (At > 1) cells we simulate this using ballistic 
transport (see the centre panel of Fig.[2]and Sec. 2.5.2[ ). 

I n regions of the grid where the cells are optically thin (At < 
1), ballistic transport becomes too diffusive and we need to resort 
to direction conserving transport or DCT (see the right panel of 
Fig.|2]and Sec. |2~5~3] l. 

We now proceed by describing these transport methods in 
more detail and show how they are combined in a general simu- 
lation. 



2.5.1. Diffuse transport 

We start with the description of conceptually the simplest form 
of transport implemented in SimpleX: at each computational cy- 
cle, which we will call swee/^jin the remainder of the text, the 
content of each nucleus is distributed equally among its neigh- 
bouring nuclei (see the left panel of Fig. [2]) We call this kind 




Fig. 2. Three principal means of transport supported in SimpleX. 
Left panel: Diffuse transport, photons from the incoming edge 
(not shown) are distributed outward along all edges (including 
the incoming edge). Centre panel: Ballistic transport, photons 
are transported along the D edges most straightforward with 
respect to the incoming direction. Right panel: Direction con- 
serving transport, photons (indicated with the dotted arrows) are 
transported as in ballistic transport but their direction is stored 
indefinitely in a global set of solid angles. 



of transport diffuse as it has no memory of direction. In a ho- 
mogeneous distribution of nuclei, the transported quantity will 
diffuse outwards, spreading spherically from the position of a 
source. This type of transport is appropriate for photons that are 
either scattered diffusely or absorbed and re-emitted in random 
directions. 



2.5.2. Ballistic transport 

Now imagine a group of photons that is transported along a 
Delaunay edge to a certain nucleus. Suppose that the nucleus 



As we will see in Sec. 



4.1 



anisotropic scattering processes can be 
simulated straightforwardly by assigning weights to the outgoing edges 
so that more radiation is transported in the appropriate directions. 

4 Although in our current implementation each photon takes exactly 
one step per sweep, this need not be the case in general. 



represents a finite optical depth. A fraction of the photons will 
be removed from the group by the interaction and another frac- 
tion will fly straight onward. Diffusive transport is not suited to 
describe this behaviour, so we introduce ballistic transport. 

In the ballistic case, the incoming direction of the photons 
is used to decide the outgoing direction (introducing a memory 
of one step into the past). In the generic Delaunay triangula- 
tion there is no outgoing edge parallel to the incoming one, so 
the outgoing photons are distributed over the D most straight- 
forward directions, where D is the dimension of the propagation 
space (see centre panel of Fig. [2]) As such, we ascertain that 
for an isotropically radiating source the complete 'sky' is filled 
with radiation because the opening angle associated with each 
edge on average corresponds to 2;t/A or 47r/A in two and three 
dimensions respectiveljP] 

Note that due to the random nature of the directions in the 
Delaunay grid, radiation will tend to lose track of its original di- 
rection after several steps, a property we will call decollimation 
as it will steadily increase the opening angle of a beam of radia- 
tion as it travels along the grid (see Fig. [3] This property renders 
ballistic transport appropriate for highly to moderately optically 
thick cells only, where just a negligible amount of radiation has 
to be transported further than a few steps. If we take unity as 
a lower limit for the optical depth of a cell for which ballistic 
transport is used, at every intersection a fraction of ( 1 - 1 /e) of 
the photons gets absorbed and the cumulative average deflection 
(decollimation) becomes 



n=\ 



(6) 



where 6d is the decollimation angle per ballistic step. In 
Sec. 3.3.1 we will measure 6 D and describe the consequences 



of decollimation in more detail. 

A visual extension of the statement given by Eq.|6]i is shown 
in Fig. [4] From the figure it is evident that the fraction of photons 
that (ever) gets a deflection larger than 45° centred around the 
initial direction falls off sharply with the optical depth of a cell. 
Only for optically thin (say 0.2) cells the fraction of photons 
whose deflection stay under 45° is lower than 0.5. 

In a realistic cosmological simulations, the effect of decolli- 
mation results in diffusion which softens shadows behind opaque 
objects (e.g. filaments and halos). The diffuse radiation field will 
penetrate into the opaque objects and ionise the high density 
gas inside. This results in ionisation of dense structures at too 
early times and the stalling of the ionisation front farther from 
the source. This clearly necessitates the introduction of a means 
of transporting photons in the optically thin regime. 

2.5.3. Direction conserving transport 

If the cells are optically thin, the loss of directionality introduced 
by ballistic transport over many steps becomes prohibitive and 
we switch to direction conserving transport (DCT for short). 
Here the radiation is confined in solid angles corresponding to 
global directions in space. So, if a photon has been emitted in 
a certain direction associated with a solid angle, it will remem- 
ber this direction while travelling along the grid. This effectively 
decouples the directionality of the radiation field from the direc- 
tions present in the grid (see right panel of Fig. 13]) 



5 Another way to look at this is that, for an isotropically emitting 
source, the number of nuclei that receive radiation must scale as the 
square (cube) of the travelled distance for two (three) dimensions. 
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Fig. 3. Example of decollimation in the plane for five bal- 
listic steps. The arrow indicates the initial influx of photons. 
According to the ballistic transport mechanism, photons are 
transported along the D most straightforward edges with respect 
to the incoming direction. The color coding is as follows: with 
every step the photons get a color that is less red and more blue. 
As the angle between adjacent edges is large in 2D (60 degrees 
on average), the bundle loses track of its original direction in 
only a few ballistic steps. 




Fig. 4. Fractions of photons (contours) with deflection angle 
within 45° of the initial direction as a function of the number 
of ballistic steps (abscissa) and optical depth (ordinate). From 
the figure it is evident that decollimation of photons is only po- 
tentially problematic when the cells are very optically thin. 



Note, however, that the transport of photons still happens 
along the three most straightforward Delaunay edges of the grid, 
where 'straightforward' is now with respect to the global direc- 
tions of the solid angle^] These straightforward edges may lead 
to nuclei that lie outside the solid angle associated with the di- 
rection of the photons. 

In this sense, we now have two types of angular resolution 
in our method. The first is related to the size of the solid angle 
in which radiation is confined spatially which is set by the num- 




6 DCT can thus be viewed as ballistic transport with complete mem- 
ory of direction. 



Fig. 5. Shadow behind a dense filament irradiated by ionising 
radiation (indicated as a black dot). The hydrogen number den- 
sity is plotted logarithmically in greyscale and ranges between 
2.8 x 1CT 5 to 0.12 cm 3 . The red and blue lines indicates the con- 
tour where the neutral fraction of hydrogen is 0.3. The thin and 
thick red lines are for ballistic and direction conserving transport 
(DCT) respectively. The blue line is obtained with the C 2 -ray 
code (Mellema et al.|2006) . The contour for DCT shows a sharp 
'shadow' (in good agreement with the C 2 -ray result) where the 
dense filament is left neutral whereas ballistic transport results in 
a more diffuse, softer, shadow. Moreover, the ballistic ionisation 
front stalls with respect to the DCT result. This is another symp- 
tom of spurious diffusion: the positive radial component of the 
diffuse radiation is smaller than it should be resulting in more 
ionisations close to the source and less flux into the ionisation 
front. 



ber of Voronoi neighbours of a typical nucleus (15.54- ■ ■ in 3D). 
We will call this spatial resolution and we stress that it depends 
solely on the nature of the Delaunay triangulation and as such is 
not adjustable. The second type of angular resolution is set by 
the global division of the sky in arbitrarily many directions (not 
necessarily a constant along the grid) and we will call this the 
directional resolution. For the tests presented here, we use 40 
directions (a directional resolution of 40) implying a solid angle 
of nj 10 sr for each unit vector. 

A technical consequence of the approach sketched above is 
the need to divide space into equal portions of solid angle whose 
normal vectors are isotropic. There are many ways to divide the 
unit sphere into equal patches but the requirement of isotropy 
cannot be met exactly in general. 

For maximal flexibility in angular resolution we have con- 
structed sets of unit vectors that are distributed isotropically in 
angle using a simulated annealing scheme. 

One major drawback of the introduction of global directions 
in the simulation is that they will give rise to artifacts much 
like the ones observed in hydrodynamical simulations on regu- 
lar grids. This in turn must be counteracted by randomisation of 
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these global directions at appropriate time intervals which makes 
DCT computationally relatively expensive. 

2.6. Combined transport 

The three means of transport described above are applied simul- 
taneously in general. If the total optical depth r tot for a given cell 
is due to multiple extinction processes 



T tot 



(7) 



the fraction, /;, of the incoming ray of photons that will be re- 
moved by process i is given by 



fi = 



T tot 



(8) 



Depending on the physical nature of this process, the photons 
are either removed from the bundle (e.g. to heat up the medium) 
or redistributed isotropically (or with some directionality) using 
diffuse transport. The remaining photons (those that are not re- 
moved from the ray) are transported with either ballistic or direc- 
tion conserving transport, depending on the total optical depth of 
the cell (ballistic if T to t > 1 and DCT if T tot < 1). The fraction of 
the incoming photons that is treated with one of the three trans- 
port methods is thus fully determined by the grid. 

3. Anisotropy and its consequences 

The 'ideal' properties of the Delaunay triangulation notwith- 
standing, the probability density distribution of the angle be- 
tween adjacent Delaunay edges is quite broad (Icke & van de 
Weygaert|[T987l |Okabe et ai]|2000] Sec.5.5.4), so that, even 
though the average triangle is the 'fattest' possible, many 'thin' 
triangles will occur. This situation may be changed by itera- 
tively adapting the underlying point process, in such a way that 
each nucleus comes to coincide with the centre-of-mass of its 
Voronoi region (see e.g. Lloyd|1982 for such an algorithm), re- 
sulting in so-called Centroidal Voronoi Tessellations (see e.g. |Du| 
et al.|1999) . This procedure produces 'most spherical' Voronoi 
regions (actually hexagons if D = 2), but is in general far too 
costly for practical computations in which the triangulation must 
be re-computed frequently. Moreover, the shifting of the nuclei 
implies that the connection with the physical properties of the 
underlying medium is no longer entirely faithful. Another rea- 
son to refrain from the use of Centroidal Voronoi Tessellations is 
the introduction of regularity and hence symmetry in the grid. In 
two dimensions, for instance, the hexagonal Voronoi cells tend 
to align, introducing globally preferential directions in the trans- 
port of radiation. This will inevitably give rise to artifacts in the 
radiation field transported on such a grid. 

Although one cannot do better than using the Voronoi- 
Delaunay construction on a random (Poisson process) point set, 
it is not obvious that this still holds when the point set is in- 
homogeneous or anisotropic. Inhomogeneity is, of course, the 
property we will encounter in all practical cases. The probabil- 
ity distribution of the directions of Delaunay edges on a Poisson 
nucleus is isotropic, but this is no longer the case when the dis- 
tribution of the nuclei is structured. 

This inherent anisotropy of the graph introduces a bias that is 
the cause of some undesirable effects. It is these effects, and their 
treatment, that we will consider here. We have tried in vain to lo- 
cate a mathematical treatment of related questions, such as: [a] Is 
the Delaunay triangulation the one that maximises the isotropy 



of the edges emanating from a given nucleus? [b] Can a process 
be found that adjusts the positions of the nuclei in such a way as 
to increase the isotropy of the edges? Concerning [a], we con- 
jecture that the answer is yes, because (at least superficially) that 
would seem to follow from the minimax property of the angles 
between the edges. As to [b], we have looked for a procedure 
analogous to the 'centre-of-mass' stratagem mentioned above, 
but have been unable to find one. 

Thus, we must face the consequences of the anisotropy bias 
on inhomogeneous point processes. We note in passing that such 
a bias will exist locally even in the case of a Poisson process for 
putting down the nuclei, because due to shot noise every instance 
of a random point process is locally anisotropic, even if it is ho- 
mogeneous and isotropic in the mean. 

3.1. Error measure 

For a homogeneous Poisson distribution of nuclei, the expec- 
tation value for the number of Delaunay neighbours (and thus 
edges), A, is 6 in two- and 15.54 • • • in three-dimensional space. 
These Delaunay edges have no preferential orientation and their 
statistical properties are well known (e.g. Okabe et al.|2 000). 

Now consider the case that the distribution of nuclei is inho- 
mogeneous. Spatial gradients then appear in the density of nu- 
clei, n(x), and the Delaunay edges connecting these nuclei are 
no longer distributed evenly over all possible orientations. At a 
given nucleus there will be, on average, more edges pointing to- 
wards high-density regions than away from them. 

This can be quantified as follows. Without loss of general- 
ity we may assume a distribution of nuclei that has a gradient 
in some fixed direction x (provided that the characteristic length 
scale of the gradient is much larger than the mean distance be- 
tween nuclei, a provision which we will assume to be fulfilled 
from now on). 

Take a cross section perpendicular to the direction of the gra- 
dient through the box at arbitrary position xq. The resulting plane 
with surface S will be pierced by Delaunay edges connecting nu- 
clei on either side of the plane. 

The number of edges piercing the plane can be estimated as 
follows. The local density of edges is the product of the number 
density of nuclei multiplied with A. An edge is able to pierce the 
plane when two requirements have been fulfilled: 

1 . Its projected length must be larger than the distance between 
its originating nucleus and the plane. 

2. It must point in the correct direction. 



The expectation value of the Delaunay edge length, £), 
bv 



by 

D = %n(x) 



■1/3 



is given 



(9) 



in three dimensions, where £ = 1 .237 ■ ■ • ( |Ritzerv eld 2007 
Eq.(3.22)). An edge emanating from a nucleus at position x 
thus will only pierce the plane if its projected length exceeds 
Ax = \x — xq\. Hence, the first requirement is fulfilled when 



Ax < D cos 6 , 



(10) 



where 8 is the angle between the edge and the direction of the 
gradient. The orientation of these edges is random^] so we must 



7 The orientation might be correlated with the direction of the gradi- 
ent but we will neglect this at the moment as it would only influence our 
results in second order. 
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average the cosine over a half sphere which yields a factor of one 
half. 

The effective length of the edges is thus D/2 which can be 
interpreted as follows. Only nuclei inside a slab of thickness £), 
centred at Xq contribute to the density of piercing edges. 

The second requirement effectively excludes half (up to first 
order) of the edges because they point away from the plane. This 
statement is equivalent to noting that every piercing edge con- 
nects exactly two nuclei at opposite sides of the plane. Including 
both factors we find that the number of piercing lines, N p , is 
given by 



N p (x) = 



n(x)ADS 



(11) 



3.2.1 . Physical slow down 

The first effect reduces the transport velocity and can be inter- 
preted as a physical phenomenon: if we identify the length of a 
Delaunay edge with the local mean free path of the transported 
quantity (e.g. photons), the shorter edges simply express the fact 
that we have entered a region of increased optical depth where 
it takes more mean free path lengths to traverse a given physical 
distance. It has been shown (Ritzerveld & Icke 2006) that identi- 
fying the average Delaunay edge length with the local mean free 
path of the relevant processes is a natural choice in the construc- 
tion of the triangulation and the observed behaviour is therefore 
both expected and physical. 



The surface density of lines piercing the slab, cr(x), is now sim- 
ply defined by 



3.2.2. Drift 



cr(x) 



N p n(x)AD 



(12) 



The sought after fractional excess, E(x), of parallel (with respect 
to the gradient) over anti-parallel Delaunay edges is thus given 
by differencing cr(x) over a sufficiently smallj interval Ax and 
division by Kn(x) 



E(x) = 



1 Acr(x) 



£ An(x) 



An(x) Ax 4n(x) 4 ^ Ax 



(13) 



where we have used Eq.|9]l to eliminate T>. The excess of edges 
pointing towards higher density region may have a signifcant ef- 
fect on quantities transported along these edges. In the next three 
sections we will describe and quantify these effects. We note in 
passing that in many practical applications (also see Sec. [6]i the 
point density n(x) of nuclei is taken to be proportional to a power 
a of the mass density p(x): 



n(x) oc p{x) a . 



In that case, Eq.( 13 1 becomes 



E(x) = 



£ A log n(x) aD A log p(x) 



4n(xy> 3 Ax 



Ax 



(14) 



(15) 



If a = 3, the length T) is proportional to the mean free path of 
the photons ( |Ritzerveld & Ick e 2006). In many respects, this is 
the 'ideal' case, because is in fact 'transport-homogeneous' as 
experienced by the photon; every step is of equal optical depth. 
However, Eq.(fT5j) shows that such an 'ideal' case has a gradient 
asymmetry that is three times worse than is the case if n(x) oc 
p(x). 

3.2. Effects on diffusive transport 



In Sec. 2.5.1 we stated that for diffusive transport in a homoge- 



neous distribution of nuclei, the radiation propagates spherically 
away from a source. On a graph corresponding to an inhomoge- 
neous distribution of nuclei, however, this is not the case. The 
spreading of the transported quantity will not be spherical any- 
more due to two effects: 

1 . The Delaunay edges are shorter when the nuclei are spaced 
more closely. 

2. The orientation of the Delaunay edges is no longer isotropic: 
more edges point towards the overdense regions. 



8 We emphasize the discrete nature of this derivative by noting that 

4n A v — fTt4ll — Cn^v^-l/^ dn 



An= fAx = D, 

dx dx 



The second effect, quantified by Eq.( 13 i, is an artifact of the 



Delaunay triangulation itself and causes unphysical behaviour. 
When there are too many edges pointing into the overdense re- 
gions, the transported particles are deflected into those regions 
and the direction of propagation tends to align with the gradient 
(see Fig.|6]for an example in the plane). We call this effect drift 
and will now proceed to quantify its consequences for diffusive 
transport. 

A photon scattering at a nucleus has the following probabili- 
ties of moving in the dense (subscript d) or underdense (subscript 
u) direction: 



Vd 

Pu 



1 E(x) 

2 + ~2~ 

1 _ E(x) 

2 ~ ~T~ 



(16) 



The expectation value do of the drift per scattering event is there- 
fore proportional to the E(x) part of an outgoing edge in the di- 
rection of the dense region. The multiplication factor, which can 




Fig. 6. Schematic example of a nucleus and its edges subject 
to a gradient in the number density of nuclei in the positive x- 
direction. More edges point toward the over-dens region (along 
the gradient). If every edge would transport an equal number 
of photons to neighbouring nuclei, the anisotropy of outgoing 
edges results in an unphysical net flow along the gradient (indi- 
cated by the arrow which is the vector sum of the edges scaled 
down by roughly a factor of three.) 



etc. 
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be found by integrating over a half-sphere, is 1 /2 since not all 
edges point exactly to the right: 



D, 



-E(x), 



(17) 



in which D is again the expectation value of the Delaunay edge 
length (see Eq.Q). 

Say we set up a scattering experiment where a number of 
photons is placed at one position of a triangulation with a den- 
sity gradient. We expect the photons to diffuse outward with an 
ever decreasing radial velocity (the distance travelled, L, scales 
with the root of the number of steps) and at the same time drift 
towards the dense region with a constant velocity. There comes a 
time (or distance) at which the drift is equal in magnitude to the 
diffusion radius. We define the drift length, Ldrift, as the scale 
at which the diffusion distance is equal to the distance trav- 
elled through drift. Roughly speaking, diffusion dominates for 
L < Ldiift and drift dominates for L > Ldrfft- Setting the diffusion 
distance equal to the distance travelled through drift, 



-E(x)N E , 



we find that the number of steps at equality is given by 
N E = 4/E 2 (x). 

Using Eq.([T3} and Eq.Q gives an equality length of 
8n(jt) 



Vn(x) 



(18) 



(19) 



(20) 



independent of D. Hence we have found a local expression for 
the relative importance of the unphysical drift. 

For a given density distribution n(x), the length L^ix) can 
be evaluated everywher^] This parameter can be interpreted as 
follows. Say the minimum L^nftix) for all x is N times the box 
side length, then the drift can be at most 1 /N of a box side while 
the radiation scatters throughout the box. In any case we must 
make sure that Ld r ift(Jt) is much bigger than the box side length. 

Thus, even if the density contrast is very small, it must not 
fluctuate too severely. Notice that it does not help to increase 
the number of nuclei, since both n(x) and Vn(x) in Eq.([20| 
scale with that number. Using more nuclei surely reduces the 
anisotropy per nucleus, Eq.( 13 I, but because the number of steps 
to be taken increases accordingly, the effect simply adds up to the 
same macroscopic behaviour. Only the shape of n{x) is relevant. 

In order to get a feel for the resulting constraint on the 
grid we consider the case for which La r if t {x) does not depend 
on position. Demanding LdnftW = constant in Eq.(20l implies 
an exponential density distribution. For such a distribution, the 
anisotropy is smeared out maximally over the domain. In this 
case, if we want the drift to be less than 1/8 of the box length 
the density contrast must be less than a factor e as 2.71 • • ■ imply- 
ing that the restrictions put by our isotropy demands are rather 
stringent. 

3.2.3. Clustering 

We have seen that the spurious drift for diffuse scattering puts 
restrictions on the density contrasts that can be simulated by the 



9 It may surprise the reader that values for E(x) and D are taken to 
be local while the argument seems to involve a domain in which these 
values could change. What actually happens is that we adopt the local 
values for the whole domain in Eq.d 1 81 



plain implementation of SimpleX introduced in Sec. ( 3.2.2 1. In all 
of this we have assumed that the scale of the density fluctuations 
is comparable to the box size. This is not actually a restriction: 
if we are interested in a case where the density fluctuations are 
of a much smaller scale than the box size, we can just put an 
imaginary box around each density fluctuation and use all the 
quantitative results from above. In doing so we see that for any 
given 'snapshot' too many photons will be present in the local 
overdense regions, but if there are no overall density contrasts on 
large scales, there will be no significant macroscopic drift. The 
effect of the local drift on small scales can still influence results 
on a large scale, however. Not only does the drift influence the 
average position of the photons, it also influences the standard 
deviation around this average. Isotropic scattering maximises the 
spreading of photons, but in an extreme case where for example 
at every nucleus 90% of the photons move in the same general 
direction, these will stick together for a longer period and there- 
fore the size of a 'light cloud' will grow more slowly. This is an 
effect that shows up if we have a highly fluctuating density field 
on small scales. A case in point is the simulation of the filaments 
of large scale cosmic structures. 

Consider a density distribution that is homogeneous in the y- 
and z -direction but highly fluctuating in the x-direction. For the 
probabilities to travel into the dense or underdense regions we 
again use Eq.([T6| but with the difference that now d and u no 
longer denote global directions. We describe the transportation 
process by a binomial distribution. If there are no large scale 
density contrasts the drift will be zero on average, but for the 
standard deviation we find: 



cr(N) = (y/Np(l-p)) 



(Vl-£ 2 «) 



(21) 



where N is the number of sweeps. For small values of E, this 
reduces to 



^(A0^(i-|(|»'(xM^ 3 | 2 )) 



(22) 



This factor is always smaller than unity, indicating a reduction 
of the spreading of the photons. This effect can be significant 
for small-scale fluctuations with either a very high amplitude or 
a very short length. The effect slowly becomes smaller if the 
number of nuclei is increased. The value of E(x) comes close 
to unity or exceeds unity only if the characteristic density fluc- 
tuation length is smaller than a Delaunay length, which we had 
pointedly excluded. 

3.3. Effects on ballistic transport 

For ballistic transport, similar problems arise as for diffusive 
transport but the fact that the transport is locally anisotropic 
complicates the picture. We first describe some properties of this 
kind of transport on a homogeneous grid in order to appreciate 
deviations from the 'normal' case later on. 

We identify two distinct phenomena that photons travel- 
ling ballistically may experience: deflection and decollimation. 
Deflection is here understood as 'loss of direction', where the 
direction is given by the vector sum of the three most straight- 
forward directions. Decollimation is defined to be the effect of 
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the increase in opening angle of a beam of photons as they are 
transported ballistically. The first phenomenon depends on the 
evolution of the vector sum of the three most straightforward di- 
rections whereas the second phenomenon is related to the angu- 
lar separation between these directions individually. In the case 
of a homogenous distribution of nuclei, the net effect of deflec- 
tion will vanish because there is no preferential direction in the 
grid. 

3.3.1. Decollimation 



As stated in sec. 2.5.2| the angular resolution, and therefore the 
minimal opening angle, is a property of the Voronoi cell and 
is thus fixed for the chosen triangulation. To exemplify this we 
consider a typical Voronoi nucleus in 3D, connected to A neigh- 
bouring nuclei. The expectation value for the solid angle, £2, sub- 
tended by each edge is thus 



Q 



An 
A~ 



71 

4* 



(23) 



In Fig.[7]the distribution of angles is given for a grid of 10 
homogeneously placed nuclei in three dimensions. For the vector 
sum the average departure from the incoming direction is about 
15° (solid line). This implies that a photon loses all knowledge 
of its original direcion in typically (90/15) 2 = 36 steps. 

For the distribution of the three separate most straightfor- 
ward edges we find a standard deviation of about 39° which 
means that after (90/39) 2 « 5 steps the photon has lost all mem- 
ory of its original direction. Note that the width of the vector sum 
of the weighted edges is smaller than that of the most straight- 
forward edge (indicated by '1' in the figure) alone. 

In general, the gradual loss of direction is not a major con- 
cern since most sources emit isotropically anyway, but the ef- 
fect becomes important when many edges are traversed, because 
then the photons behave diffusively, producing an intensity pro- 
file that deviates from the correct r~ (flL1) -form. Simulations in 
which the mean free path for scattering or absorption is much 
smaller than 5 edges are fine in this respect, because then the 
ballistic photons never enter this random walk regime (see also 
the discussion in Sec. [7]). To observe the implications of this sta- 
tistical statement in the transport of photons we construct a trian- 
gulation of a random point distribution containing 10 5 homoge- 
neously placed nuclei in a cube of unit dimension. In the centre 
we define a small spherical volume with radius 0.05 and place 
a number of photons in each nucleus contained in this volume. 
The photons are assigned a direction by sending them along the 
edge whose direction is maximally parallel to one of the coordi- 
nate axes (we will choose the positive x-axis). So, when applying 
our ballistic transport method, the photons will start to move in 
the positive x-direction with a step-size equal (on average) to 
cos 6d where Op is the mean decollimation angle. Every step 
means a multiplication by a factor cos do, until after infinitely 
many steps the effective step size in the x-direction is zero. In 
Fig. [8] the step-size as a function of sweeps is shown together 
with a fitting function of the form 



f(N s )= A cos(0 D ) 



N> 



(24) 



where N s the number of sweeps. We find that in the case with- 
out ballistic weights (filled squares; see Sec. |4.2| for the descrip- 
tion of ballistic weights) G D = 0.67 ± 1% (39° + 1%) in accor- 
dance with the value obtained from the grid statistics presented 
in Fig. [7] 



3.3.2. Deflection 

We are now ready to quantify the effect of anisotropy of the trian- 
gulation on the ballistic transport of photons, which travel along 
the three edges closest (in angular sense) to the incoming direc- 
tion. For the sake of simplicity we will estimate the deflection 
for radiation travelling along only one most straightforward edge 
and discuss the applicability to three edges afterwards. 

Consider photons streaming perpendicular to the gradient di- 
rection (see Fig.[9]for the geometry of this situation.) The stan- 
dard deviation of the deflection of the outgoing edge with respect 
to the incoming direction is typically 15° (see Fig. [7} for a ho- 
mogeneous grid but will in general depend on the local value for 
the gradient. Referring to Eq.(23l, the expected opening angle 
depends on the effective number of outgoing edges in that di- 
rection. In other words, as the anisotropy increases the number 
of edges pointing toward the over-dense region, the angular res- 
olution increases accordingly. This motivates one to define the 
direction dependent number of outgoing edges, A e ff , as 



A eff (0) = A[l + E(*)cos 0], 



(25) 



where tp is the angle with the gradient direction. If we interpret 
the solid angle as a projected circular area on the unit sphere 
we can estimate the maximal deflection angle, for ballistic 
transport as 



\Aeff(0)' 



(26) 



For a homogeneous distribution of nuclei, this angle evaluates to 
approximately 30° with a typical value of 15° as expected from 
the analysis shown in Fig. [7] If we include this empirical factor of 
one half in Eq.(26 1 and approximate the arcsine by its argument 



(correct within 1% for angles smaller than tt/12) we obtain 



(A[l +E(x)cos<p]) 



-1/2 



(27) 
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Fig. 7. Normalised distribution of angles between the incoming 
direction and the vector sum of the (three) most straightforward 
pointing edges (solid line) and the separate most straightforward 
edges (dotted, dashed and dot-dashed lines) for a homogeneous 
distribution of nuclei. The distribution of the most straightfor- 
ward edges added is shown as the long dot-dashed line with label 
'Added' and is relate d to the decollimation effect. The standard 
deviation (FWHM/ Vln256) is about 39° for the added edges 
and 15° for the vector sum of the edges. 
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Fig. 8. The projected step-size in the x-direction as a function 
of sweeps together with a least squares fit for the unweighted 
ballistic transport (solid line) and with ballistic weights as de- 
scribed in Sec. |4.2| (dashed line). A region indicating unit stan- 
dard deviation is shaded in grey around the markers. Error bars 
corresponding to cr/ VA - 1 where cr the standard deviation and 
N the number of experiments executed. For the uncorrected case 
the fitted value for the deflection angle do in Eq. ( 24 1 is 39° ± 
1%. For the weighted case we find 9d = 26° + 3%. 
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Fig. 9. Geometry of radiation travelling perpendicular to the gra- 
dient direction. The deflection angle toward the under-dense re- 
gion, 6 U , scales with the radius of the area S u and similarly for 
the deflection angle towards the over-dense region, 9^. 



where we have used Eq.( 25 1 to substitute for A e ff . This result can 
be approximated to first order by 



a -i/2m _ E(x)CQS( 



■]■ 



(28) 



To obtain the deflection per step over the grid, we average 0^ 
over azimuthal angle while projecting along the direction of the 
gradient 



1 C\ 

2?r Jo 



6& cos 



E(x) 
4VX' 



(29) 
(30) 



The sign of 6 e ff is negative, which means that the resulting de- 
flection is into the region of lesser density. This result might sur- 
prise the reader at first. One could naively predict the photons to 
deflect into the direction of higher density (similar to the diffuse 
drift) but the situation is exactly reversed for ballistic deflection. 

Note that this effective deflection angle is an upper limit for 
the deflection encountered in a simulation because of the follow- 
ing. The effect is maximal for radiation travelling perpendicular 



to the direction of the gradient and this is the situation we have 
used as a starting point for the above derivation of 6 e s. Secondly, 
the fact that 6*d < U will favour the selection of edges in the over- 
dens region effectively diminishing the deflection, an effect that 
we have neglected in the derivation above. To include this effect 
one would have to know the distribution of outgoing edges as a 
function of angle with the gradient, and assign probability to the 
selection of an edge accordingly. As we shall see in Sec. 5.3 the 
omission of this effect does not seem to be of much importance 
for our predictions. 

Sending radiation along three edges rather than one will fur- 
ther decrease the effect of deflection. This is expected as the de- 
flection of one edge is larger than that of the vector sum of three 
outgoing edges as we saw in the case of a homogeneous grid. 

The growth of the deflection when traversing the grid can be 
found by differencing Eq.(30 1 by the typical step-size, T> 

Afleff = 1 V»(X) 

A/ 16 VX n(x) ' 



Analogous to the drift length of Eq.(20 1 we can define a de- 
flection length, Ldef , at which the cumulative deflection is equal 
to say, 7r/4, 



^def - 



47r VAn(jc) 
Vh(jc) 



50 



n(x) 
Vn(x) ' 



(32) 



Comparing with Eq.(20l we see that both the diffuse drift and 
the ballistic deflection give approximately the same restrictions 
on the box size, the only difference being a factor 50 instead of 
8. Again the interpretation is simple. If Ldef(*) is strictly larger 
than say five times the box size, the deflection will be less than 
one-fifth of tt/4 in travelling across the box in a 'straight' line. 
Remember that this effect is added to the expected decollimation 
from Sec l3~3~Tl 

4. Weighting schemes 

In this section we describe several possibilities to correct for the 
unphysical effects described in Sec. [3] 

4. 1 . Diffuse transport 



A straightforward cure for the problems addressed in Sec. 3.2.2 
and Sec. 3.2.3 (diffuse drift and clustering respectively) is to as- 
sign weights Wj to the edges emanating from a given nucleus in 
such a way that the anisotropy vanishes. This means that the frac- 
tions of the quantity transported to the neighbours are no longer 
equal to l/N but directly proportional to the solid angle that the 
corresponding Voronoi face spans. We refer to Fig. 10 for an ex- 
ample in the plane, the three-dimensional case is analogous. We 
have explored this possibility and implemented three different 
weighting schemes in our method, each with its specific virtues 
and drawbacks: 

1 . Voronoi weights: based on the natural properties of the trian- 
gulation. Advantages: automatically and 'naturally' adapted 
to the physics of the transport problem. Disadvantages: in- 
herently noisy. 

2. Icosahedron weights: based on a division of the unit sphere 
using the icosahedron. Advantages: flexible (does not de- 
pend on type of triangulation), refinable. Disadvantages: 
computationally expensive. 

3. Distance weights: based on the distance to a neighbour 
squared. Advantages: very fast. Disadvantages: Only empir- 
ical evidence of its correctness. 
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Fig. 10. Solid angles for two edges in the case of the Voronoi 
weighting scheme (cell 3) and the icosahedron weighting 
scheme (cell 1). In the case of the Voronoi scheme, the weight 
assigned to a Delaunay edge corresponds directly to the magni- 
tude of the solid angle of its corresponding Voronoi cell. In the 
icosahedron scheme, the weight is proportional to the solid an- 
gle a Delaunay line occupies considering the angular vicinity of 
its neighbouring edges (in 2D the solid angle is thus bound by 
the two bisectors shown as dotted lines in the figure.) 




Fig. 11. Geometry of the COG weighting procedure. In general, 
the vector sum of all (weighted) Delaunay edges (black arrows) 
is non-zero (red arrow). The Delaunay edges j and i are most 
parallel and anti-parallel to this vector sum respectively. Their 
weights are adjusted such that the vector sum vanishes in its cur- 
rent direction. To this end, the vector sum and its reflected coun- 
terpart (green arrow) are projected on j and i, indicated with the 
smaller black arrows. Half of this projection is added to edge i 
and half is subtracted from j. The resulting vector is shown in 
blue and it exactly cancels the vector sum. The adjustment of 
weights can still result in a non-zero (although smaller) vector 
sum in some other direction. The above procedure is repeated 
iteratively until the norm of the vector sum is smaller than a pre- 
defined tolerance. 



4.1 .1 . Centre of gravity weighting 

In all three weighting methods, the weights sum up to unity, 
guaranteeing conservation of photons. In addition the vector 
sum of the weighted Delaunay edges emanating from a nucleus 
should be zero to conserve momentum of the radiation field. We 
obtain this by adjusting the weights of the (two) edges most par- 
allel and anti-parallel to the 'centre of gravity', of the nucleus 
defined by 



z 



win 



(33) 



such that in that direction the magnitude of *P vanishes (see 
Fig. 11) Here r, denotes the direction of the i-th Delaunay edge. 



We repeat this procedure until the absolute value of *P is smaller 
than some given tolerance. In our experience, this procedure 
reaches a (relative) tolerance of 10~ 7 within twenty iterations 
for the Voronoi method and about ten for the Icosahedron and 
Distance weighting schemes. 

4.1.2. Voronoi weighting 

The Voronoi weighting scheme uses the faces of the Voronoi 
cells to calculate the w,- for each nucleus by giving an estimate 
for the solid angles subtended by the walls of the Voronoi re- 
gions, normalised to unity. The weights are then 



Wi = 



Sid; 



2 ' 



(34) 



where S , denotes the surface of the Voronoi face perpendicular 
to the Delaunay edge connecting the current nucleus and its 2-th 
neighbour and dj is the length of that edge (because the Voronoi 



face cuts the Delaunay edge in the middle, 0.5ii, is the distance 
from the nucleus to the face). Note that Eq.([34| is an approx- 
imate expression neglecting projection effects for large angles. 
With some computational effort, this estimate could be refined. 
The solid angle subtended by a Voronoi face depends on the size 
of the cell but also on the precise position of the nuclei. Referring 
to Fig. 



10 



we see that cell '3' has a relatively small surface be- 
cause its nucleus is close to that of neighbour '4' and its nucleus 
is a bit further from the central nucleus than that of neighbour 
'4'. If that last statement were reversed, neighbour '4' would 
have the smaller surface and would thus get the smaller weight. 
This property may seem harmful at first but it is of stochastic 
nature, meaning that some noise will be introduced but no sys- 
tematic error. 

To demonstrate our method we have constructed a Delaunay 
triangulation of 10 5 nuclei in three-space with a linear gradient 
in n{x) along the horizontal direction which runs from 0.005 on 
the left to 0.995 on the right. As a measure for the anisotropy 
in orientations of Delaunay edges, we take the angle between an 
edge and the direction of the gradient and plot the number of 



edges in an angular bin (see Fig. 12 1. Because we are interested 



in the relative deviation from a horizontal line (which would in- 
dicate the isotropic situation), the results are given as a fraction 
of the average, f(6). The data displayed in the topmost panel 
corresponds to the weighting schemes described above which 
can be thought of as a basic correction plus a refinement thereof 
(the 'centre of gravity' correction). We show the effects on the 
angular distribution of the edges as we apply these corrections 
cumulatively. 

The initial anisotropy apparent from the inclination of the 
line labeled 'Uncorrected' is reduced significantly (to about half 
its original value) by application of the Voronoi weighting alone 
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(dashed line). After application of the 'centre of gravity' cor- 
rection, the scatter around the isotropic value of unity is below 
0.25% except at the outer edges where the normalisation blows 
up the errors. 

Apart from the anisotropy due to the gradient, the triangula- 
tion shows some noise of order 0.5% itself (the light grey area 
around the uncorrected lines in the three panels have a width of 
one standard deviation). Moreover, after all weighting has been 
applied, some noise remains (a unit-standard deviation area is 
shaded in dark grey). This noise is related to features of the tri- 
angulation itself. If the triangulation itself has more edges in a 
certain direction, this cannot be corrected with a local weighting 
scheme as it is a global property subject to chance. Such noise 
can be reduced by either placing more points, or, equivalently, 
constructing several instances of the same triangulation and av- 
eraging over the results. 

Because all the information needed for the Voronoi weight- 
ing scheme is inherent in the tessellation and its triangula- 
tion, the computational overhead associated with it is potentially 
small. Unfortunately, in most tessellation software, the areas of 
the Voronoi walls are not computed with the other properties 
of the tessellation. Calculating these areas is costly as in three 
dimensions the walls are generally irregular polygons with M 
vertices where M > 3. 

4.1.3. Icosahedron weighting 

The second method is based on an 'independent' division of 
the unit-sphere into M (approximately) equal parts. We take the 
M vectors (originating from the nucleus under consideration) 
that point to these parts and assign weights w,- to the outgoing 
Delaunay edges as follows. We take a vector and calculate the 
N dot products with the Delaunay edges. The Delaunay edge 
that has the smallest dot product gets a fraction 1 /M added to its 
weight. 

We have chosen the icosahedron as the basis for our weight- 
ing scheme. We use vectors to the middle of its 30 edges together 
with its 12 vertices, yielding 42 reference vectors. 

One significant advantage of the icosahedron method is that 
it is independent of the nature of the triangulation and is there- 
fore very flexible. It becomes more accurate as the division 
of space is refined (M is increased) by taking other tessella- 
tions of the unit-sphere. The computational cost, however, scales 
linearly with M forcing us to trade off accuracy with speed. 
Unfortunately this 'independence' has a drawback as well: the 
reference vectors are oriented statically in space, introducing a 
systematic bias in those directions. In order to reduce this effect, 
one is forced to add some random noise to the procedure, for ex- 
ample by applying random rotations of the whole icosahedron, 
another costly operation. 

After the correction with icosahedron weights (see middle 
panel of Fig. 12 dashed line), the anisotropy has decreased be- 
low the 0.5% level. This immediately shows the strength of this 
method over the more noisy Voronoi scheme (also compare the 
unit standard deviation regions in dark grey). 



4.1.4. Distance weighting 

Empirically we have found that in 3D the square of the distance 
between nuclei is also a robust estimator of w,-. This quantity is 
readily calculated and is by far the fastest solution to the problem 
at hand. After the initial correction (see bottom panel of Fig. 12 
dashed line) the anisotropy of edges has been diminished to val- 
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Fig. 12. Fraction of edges, /(#), with angle n — 6 with re- 
spect to the direction of the linear gradient as a function of 9 
for different cumulatively applied corrections with the Voronoi 
(top panel), Icosahedron (middle panel) and Distance weighting 
scheme (bottom panel). The results are averaged over 100 differ- 
ent realisations of the grid to suppress shot noise. Unit standard 
deviation regions around the 'Uncorrected' curve and the final 
result are shaded in light and dark grey respectively. As a result 
of the anisotropy of the triangulation more edges point towards 
the overdense region (to the right) in the uncorrected case. Note 
that due to normalization, the results at the extreme ends are sub- 
ject to noise. 



ues lower than 0.5% even slightly better than in the Icosahedron 
scheme. The deviations from the ideal isotropic case after the 
COG correction are of the same order as in the Voronoi and the 
Icosahedron case. 

We have not found a valid explanation for the success of this 
method. Intuitively, one would think that the opening angle, Q, 
of a Voronoi wall with respect to its nucleus would scale as r~ 2 
rather than r 2 . In order to find a mathematical reason for the 
proportionality between Q and r 2 one is forced to delve deep 
into the field of computational geometry, an endeavour clearly 
beyond the scope of this text. 
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4.2. Ballistic weighting 

Assigning weights to the outgoing edges as described above re- 
moves the drift and clustering problems for diffusive transport 
described in Secf5] 

The effects of anisotropy in the ballistic case are intrinsically 
more challenging to correct because it is not a priori clear how 
the anisotropy of all outgoing edges of a nuclues should be used 
to choose three weights for the outgoing edges. 

As described in Sec |3.3.ll and Sec |3.3.2 in the ballistic case 
we must distinguish between decollimation and deflection of 
which the first one dominates the overall 'loss of direction' . 

In order to diminish the decollimation of a beam, we can 
assign weights to the most straightforward pointing edges. If we 
make these weights somehow proportional to the inner product 
of the edge and the incoming direction, the 'straightest' edges 
get most of the photons. 

There are in fact many possibilities to assign weights to the 
D most straightforward edges, optimising different aspects of the 
transportation process. If the exact direction is important, for in- 
stance, the weights should be chosen such that the vector sum 
of the resulting edges points straight ahead. On the other hand, 
minimising the decollimation will in general yield a different set 
of weights, those that maximise the length of the vector sum. 

From this myriad of possibilities we choose a very simple 
and computationally cheap approach. Edge / gets weight Wj 
given by 



Vector sum 
Added 



cos(6j)/sm(6j) 
Zjcos(0j)/sm(0j)' 



(35) 



where Oj is the angle between the incoming direction and edge 
j. Division by the sine of the angle places extra emphasis on the 
edges with smallest inner product. 

Application of this weighting scheme on the point distribu- 
tion used to produce Fig.|7]yields the distribution of edges shown 
in Fig. 13 The ballistic weights as specified by Eq.([35| have the 



desired effect of decreasing the width of the distributions of the 
edges used in ballistic transport. 

In Fig. Revolution of the step-size of ballistic transport in the 
weighted case is shown with crosses. The decollimation angle 
per step is substantially smaller than in the case without weight- 
ing (26° and 39° respectively). This decrease in decollimation 
angle consequently relaxes the requirement that the number of 
ballistic steps must stay below 5. Including ballistic weights typ- 
ically allows for up to 12 ballistic steps before the original direc- 
tion has been lost. 

Lastly, there are situations where the anisotropy of edges em- 
anating from a cell cannot be solved by any weighting scheme. 
If there is no edge in a given direction, the radiation cannot go 
there, whatever the weights are. Loosely speaking, one should 
start off with a reasonably isotropic triangulation in order to ap- 
ply a weighting scheme in a useful way. In the case of Delaunay 
triangulations in three-space, this is almost always the case. 

4.3. OCT throughout the grid 

A conceptually different solution for the problems described 
in Section [3] is to use the direction conserving transport from 
Sec. 233] throughout the simulation domain. As the angular di- 



rection of the radiation is effectively decoupled from the grid in 
this approach, drift, clustering, decollimation and deflection are 
solved for at once. The cone of nuclei that will receive radia- 
tion will however still be dependent on the grid itself as already 
pointed out in Sec. [233] 
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Fig. 13. Normalised distribution of angles between the incom- 
ing direction and the vector sum of the (three) most straightfor- 
ward pointing edges (solid line) and the separate most straight- 
forward edges (dotted, dashed and dot-dashed lines) for a homo- 
geneous distribution of nuclei with ballistic weights. The distri- 
bution of the most straightforward edges added is shown as the 
long dot-dashed line with label Added' . The standard deviation 
(FWHM/ VhT256) is about 26° for the added edges and 16° for 
the vector sum of the edges. Note that the width of the distri- 
bution of the single most straightforward edge poses a definite 
lower limit for the width of the Added' lines. 



5. Numerical examples 

Having considered the various systematic effects inherent to the 
SimpleX algorithm, we now present a number of sample cases 
demonstrating them. In general, all of the effects described above 
occur simultaneously, but for transparency we analyse them in 
isolation. 



5.1. Effects on diffusive transport: drift 

We now proceed by describing the numerical experiment de- 
signed to show the effects of diffuse drift. The simulation domain 
of unity volume in three dimensions is filled with 5 10 5 nuclei 
subject to a gradient in the point density of the form n(x) = x. 
We choose this linear form because it locally approximates every 
other type of gradient. The number of nuclei results in roughly 
32 steps across the box along the gradient direction and 38 in the 
directions perpendicular to the gradient. These numbers allow 
for most of the photons to travel through the box for more than 
100 steps before being captured in the absorbing boundaries. 

At the site closest to the centre of the domain, a number of 
photons is placed. Neither the outcome nor the speed of our sim- 
ulation depends on this number as we use floating point numbers 
to represent phtons. 

Photons are transported over this grid using diffusive trans- 
port without absorption. With every sweep the photon cloud is 
expected to grow in size. In the case of a homogeneous point 
density the photons will be distributed normally, as must be ex- 
pected for pure diffusion. The gradient in the point density will 
distort the form of the distribution function because of two rea- 
sons. First, the mean free path on the grid, D, scales with the 
point density according to Eq.(j9j alowing photons to diffuse 
faster into the underdense regions where the step sizes are larger. 
Second, the drift phenomenon described in Sec. |3.2.2"] will coun- 
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teract this physical diffusion and will move the cloud into the 
overdense region. 

In order to separate these two effects, we have performed a 
one dimensional Monte Carlo simulation of 2 10 6 'random walk- 
ers' that take steps with a size given by the recipe of Eq.(|9]i. 
Essentially, the experiment described above is emulated in one 
dimension and without the Delaunay grid as underlying struc- 
ture. The random walkers are thus expected to experience the 
physical diffusion into the underdense region only. The unphys- 
ical drift is exclusively due to the Delaunay grid and will not 
show up in the results. 

In Fig 14 the intensity weighted position of the photon cloud 



(along the direction of the gradient) versus the number of sweeps 
in the simulation is shown for SimpleX and the Monte Carlo ex- 
periments both in the case with and without weight^] We can 
see that the behaviour conforms to our expectation. The weight- 
ing corrected SimpleX result coincides with the Monte Carlo as 
shown in the bottom panel. 



Furthermore, we have taken our drift description Eq.( 13 i and 



applied it to the Monte Carlo experiment thus introducing a drift 
toward the overdense region as would be experienced by a pho- 
ton in SimpleX. The effect counteracts the physical diffusion into 
the underdense region resulting in a positive slope for the posi- 
tion of the photon cloud as a function of sweeps (see top panel of 
Fig 14 ) This procedure thus provides a direct quantitative check 



to the correctness of Eq.( 13 1 



The slopes of the Monte Carlo experiment agree within a 
thousandth of a degree with the SimpleX results both for the un- 
corrected and weighted results. This tells us several things: First 
that the weighting schemes discussed in Sec.|4]do not only cor- 
rect the orientation of the edges in a statistical sense (as shown 
in Fig. 12 1 but also allow for correct transport over these edges. 



This may seem a trivial statement but it could well be that al- 
though everything seems fine in a global sense, local anomalies 
would prevail. The transport of photons however, depends vitally 
on the local correctness of the weighting scheme and is thus a 
more stringent test. Furthermore, the recovery of the SimpleX re- 
sults with the Monte Carlo experiments suggests that Eq.( 13 1 de- 



scribes the drift phenomenon accurately in a quantitative sense. 



5.2. Effects on Diffuse Transport: Clustering 

As seen in Sec. 3.2.3| the expansion of a photon cloud is stalled 
by small (relative to the simulation domain) scale gradients in 
the grid. To show this effect we have set up the following ex- 
periment. The simulation domain is again given by the unit cube 
with a large number of photons in the centre. The 5 10 5 nuclei are 
distributed according to the following probablity distribution: 



n(x) = 0.2 + 0.8 sm\co\x\) 



(36) 



where u = 20. The density of nuclei thus inhibits concentric 
variations with a amplitude of 0.8 on a homogeneous back- 
ground of 0.2. 

We can now prepare a similar simulation with a homoge- 
neous distribution of nuclei where the box-size measured in units 
of T> is the same. Comparison of the spread of the photon cloud 
in the homogenous setup with the one according to Eq.( 36 1 will 
be a direct measure of the effect described by Eq.(22i. 

Application of a weighting scheme for diffusive transport as 
described in Sec. |4] should remove the difference between the 



10 We have used Icosahedron weights in this case but this choice does 
not influence the results. 
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Fig. 14. Intensity weighted position of a photon cloud released 
in the centre of a simulation domain with a linear gradient in the 
number denisty of nuclei as a function of the number of sweeps. 
Results obtained with SimpleX are shown with filled symbols 
and results obtained by the Monte Carlo experiments are indi- 
cated with open symbols. We use triangles and squares in the 
corrected and uncorrected case respectively. 



spread of the cloud in the two cases described above. As a check 
we will also compare the results to the analytical result given 
and the expected spread of the photon cloud in the 



by Eq.(21 



homogenous case (which is simply T) yN /2) 
In Fig. 



15 the spread of the photon cloud in terms of nor- 



malised standard deviation of the intensity weighted positions 
is shown for the inhomogeneous point density given by Eq.([36| 
(solid lines) and the analytical prediction (dashed line). The data 
are obtained from twenty runs with different realisations of the 
grid. 
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Fig. 15. Normalised standard deviation of the intensity weighted 
positions of a photon cloud expanding in a distribution of nuclei 
with small scale gradients with and without weights (thick and 
thin solid lines respectively) and the analytical expectations ac- 
cording to Eq.(21 1 (dashed line). A dotted line at cr/cr hom = 1 is 



included to guide the eye. 



As expected, the photon cloud expands (almost 6%) slower 
when no weights are used in the case with small scale gradients 
(thin solid line). Application of the weighting scheme results in 
a cloud that (after a slow start) conforms to the expected size 
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(dotted line at cr/cr hom = 1). The analytical result predicts this 
behaviour within 3% after 30 sweeps and increasingly worse for 
less sweeps (although the deviation stays under 10%). The larger 
deviations for fewer sweeps are due to the fact that the size of the 
cloud (expressed in standard deviations) is small and the results 
are divided by this small number. Note that the unit standard de- 
viation region for the weighted case is less wide for the weighted 
case (light gray) than for the unweighted case (dark gray). This is 
expected as the weigting scheme corrects for local anisotropies 
and thus reduces noise in the diffusive transport. 

We have verified that for the homogeneous case application 
of a weighting scheme does not alter the expansion speed of a 
photon cloud significantly, as expected. We can thus conclude 
that, although clustering can impose a substantial effect on the 
expansion of diffuse radiation, our weighting scheme appropri- 
ately corrects for it. 

5.3. Effects on ballistic transport: deflection 

In order to study if anisotropy has any effect on the direction of a 
beam of photons, we set up a triangulation of a random point dis- 
tribution containing 5 10 5 nuclei placed in a cube of unit dimen- 
sion with a linear gradient along the x-direction. In the centre we 
define a spherical volume with radius 0.05 and place a number 
of photons in each nucleus contained in this volume. The pho- 
tons are assigned a direction by sending them along the edge 
whose direction is maximally perpendicular the x-axis. The sim- 
ulation has been executed with ballistic transport using both one 
and three outgoing edges and additionally with the direction con- 



serving implementation of SimpleX (as described in Sec. 2.5.3 1. 
Every run has been repeated 10 times with different instances of 
the grid to suppress shot noise and obtain error estimates. 
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Fig. 16. Mean angular directionality of a photon cloud as a func- 
tion of the number of sweeps. In this simulation a linear gradient 
in the density of nuclei is present and the angle is measured rel- 
ative to the direction of the gradient, so smaller values for the 
angle are pointed toward the denser region. As expected, the di- 
rection conserving implementation of SimpleX (as described in 
Sec. 2.5.3| indicated by filled squares) has no long-term loss of 
directionality although it shows some oscillatory behaviour that 
dampens with time. Ballistic transport with one straightforward 
edge (open triangles) shows a deflection towards larger angle 
with number of sweeps of order 5 10~ 4 n rad/sweep in accordance 
with predictions. Using three outgoing edges (filled triangles) di- 
minishes the effective deflection by a factor of five. Unit standard 
deviation regions are shaded in grey for the three simulations. 



When applying our ballistic transport method, the photons 
will undergo small angular deflections into the under-dense re- 
gion as predicted in sec. 3.3.2 The expected effective deflec- 
tion for this setup is, using Eq.(30i, # e ff = 510~ 4 7r rad/sweep. 
Thus yielding a cumulative result of approximately 0.025 n rad 
after 50 sweeps. This estimate is nicely consistent with the re- 
sult from the simulation with one straightforward edge denoted 
by the open triangles in Fig. 16 As expected, sending radiation 
along three edges decreases the deflection and fluctuation of the 
angular direction in general because radiation is distributed over 
many (3"») edges, quickly reducing shot noise. When using the 
direction conserving implementation of SimpleX (filled squares), 
the radiation does not suffer from ballistic deflection. 



6. Application 

In this section we will show how to apply the results from the 
analysis above to a more (astro-)physically relevant problem. In 
short, we ask ourselves the question: given a density (or opacity) 
distribution, how do we construct a Voronoi-Delaunay structure 
that keeps the described problems acceptably small while repre- 
senting the physical problem as accurately as possible. This will 
often mean maximising the dynamic range with the prerequisite 
that low density regions still have enough resolution. 

6.1. Sampling function 

The trade-off is thus between maximising dynamic range and 
minimising the strong gradients this inevitably introduces. In 
other words: we must investigate the effects that different sam- 
pling functions have on the triangulation and its transport prop- 
erties. These functions translate the physical medium (expressed 
as an opacity or density map) to a point distribution n{x) whose 
points form the nuclei of the triangulation. As mentioned in 
Sec |3.2.2| the average Delaunay edge length can be associated 
with the mean free path of the photons in a very natural way. 
This is achieved when n(x) scales with the opacity (or density) 
to the D-th power ( [Ritzerveld & Icke|2006| Eq.( 14)) 



n(x) = <&*p D (x) 



(37) 



where <L> denotes a homogeneous Poisson process, D is the di- 
mension of the propagation space, and * denotes the convolution. 



One drawback of using Eq.(37i is that the resulting point 



process, given a substantial density contrast, places many nuclei 
in dense regions and very few in underdense regions. For many 
applications, simulating the true mean free path may thus lead to 
a required number of nuclei that is prohibitively large. 

We therefore adopt a more flexible sampling function which 
behaves differently for the extremal densities of the grid. 
According to this function, the number density of nuclei that 
generate the Voronoi grid is given by 



n(x) = 2d) 



(y- a +y 



(38) 



where y = p/po and po is a reference density that marks the 
transition between the two regimes of sampling power. In our 
applications, a will be always smaller than D (and even smaller 
than unity) in order to prohibit the over-emphasis on high density 
regions. 

According to Eq.i 



20 



i and Eq.p2|), the quantity that links the 
properties of the grid (density of nuclei and gradients therein) 
to the systematic effects described in the previous sections is 
Q n = n/Vn. We are now in the position to define, analogous 



16 



C. J. H. Kruip, J. -P. Paardekooper B. J. F. Clauwens, & V. Icke: Mathematical properties of the SimpleX algorithm 



to Q„, the quantity Q y = y/Vy which can be measured over 
the physical density (opacity) field. As the sampling function of 
Eq.(|38]l translates the physical field to the density of nuclei, the 
measured value of Q y and the upper limit of Q n (posed by the 
maximally acceptable value for either L eq or Ldef) will constrain 
the sampling parameters a and po. 

6.2. Cosmological density field 

Let us now apply this idea to a more realistic example and see 
what this all means in practice. A typical application wherein 
the strength of the SimpleX algorithm can be fully deployed is 
the epoch of re-ionisation (see Paardekooper et al. (2009 1 for 
an account of relevant test cases). Both the wide range in den- 
sities (four to five orders of magnitude), and the large number 
of sources (hundreds to thousands in realistic cosmological vol- 
umes) are treated naturally by SimpleX due to the adaptive nature 
of the Delaunay grid and the property that computational effort 
does not increase with increasing number of sources. 

We will now take a relatively small (500/h Mpc comoving) 
cosmological volume (as used in Test 4 of |Iliev et al.| f2006)) and 
construct a representative point set using constraints from anal- 
ysis of the original data. In this specific case the cosmological 
density field was presented as a regular grid of 128 3 equal size 
cells. 



6.3. Constraints 

Given the regular grid, the quantity Q y can be determined for 
every pair of adjacent cells. The thus obtained mean value for 
£ y is 1.6. 

This value implies that if we would use a linear translation 
between the mass density of the original data and the number 
density of the nuclei of the Delaunay grid, the equality length 
(see Eq. ( 20 1) would become 8 ■ 1 .6 = 13 and the deflection length 
of Eq.(|32[i becomes 50 ■ 1.6 = 81. This means that radiation can 
travel 13 box-lengths before diffuse drift starts to dominate the 
transport of radiation, and 8 1 box-lengths before ballistic deflec- 
tion reaches a cumulative magnitude of 7r/4. 

Such a linear sampling could be performed using a — 1 in 



Eq.(14i or, equivalently, a — 1 and lim Po ^o in Eq.(38i. Linear 



sampling is hardly ever desirable in a cosmological setting, how- 
ever, as it tends to place the bulk of the nuclei in high density 
regions resulting in over-sampling (many nuclei per grid cell of 
the original grid) of filaments and clumps and under-sampling of 
voids. 

Let us now proceed by finding optimal values for our sam- 
ple parameters, given the measured value of Q v and a chosen 
value for the maximal allowed equality length. An upper limit 
for the equality length, through Eq.(20i, immediately implies a 
maximal Q„. 

If, for example, we require L eq = 40 (and consequently 
Ldef = 250) we see that we need to choose a and po such that 
Q n = 5. From Eq.( 38 1 we find the expression that links the mea- 



sured value for Q Y to Q„ and thus to the sample parameters 



Qn = Qy 



y +y 



y{ay- a - x + Dy~ D ~ l ) 



(39) 



Solutions for Q n - Q y exist only for certain pairs of a and po as 
can be seen from the limit of Eq.( 39 1 



So, if we would choose a at its maximal value of Q y /Q 1 ,'"', we 
would effectively force po — > 0, resulting in the limit where all 
density regimes are sampled with the same power, a. Choosing 
a to be smaller but close to its extremal value maximises the 
dynamic range in the high density regions and allows po to at- 
tain values where still a substantial volume of the density field 
(namely the voids) is sampled with the advantageous power D. 
Solving Eq.(40 1 for Q y - 1.6 we find a = 0.32 as the maxi- 



mal value. Taking the somewhat smaller value of 0.3, the cor- 
responding value for po is 3.7 x 10~ 5 ctrT 3 (well larger than 
the smallest densities in the cosmological field which are about 
2.4 x lO' 5 ). 



6.4. The grid 



Q'r = lim Q„ = ^ 

;y->°° a 



(40) 



In the right panel of Fig. 17 a cut at z = 0.5 through the sam- 
pling obtained with these parameters is shown. The middle panel 
shows the point distribution one gets when sampling with a — 3 
throughout the grid. Evidently, this distribution does little jus- 
tice to the low-density regions while placing almost all points in 
the dense filaments and their intersections whereas the 'hybrid' 
sampling method places sufficient resolution in both the low- 
and high-density regions. Moreover, for the a = 3 sampling, 
the cell size in the filaments and dense clumps is many times 
smaller than the that of the original regular grid. No information 
is carried by this extra resolution, however, effectively wasting 
computational resources. As an extra constraint (apart from L eq 
and Ldef) one could demand that the smallest Voronoi volume is 
equal to the resolution of the original data. This could either be 
the highest AMR resolution or, in the case of SPH, the minimal 
smoothing length. 

6.5. Speed considerations 

Although we will not delve into the details here, the speed of 
SimpleX also depends on the form of the grid. Because in every 
sweep all nuclei are visited exactly once, the execution time of 
the method scales linearly with the number of nuclei. In turn, the 
number of sweeps needed to obtain a converged roughly scales 
as the number of sweeps it takes to enable every point within the 
simulation domain to 'communicate' its physical properties with 
every other other location by exchange of photons. 

Now consider the point distribution in the middle panel of 
Fig. [17] The resulting grid will have many short edges in the high 
density clump in the east of the domain. For diffusive transport or 
ballistic transport over many (more than roughly 5) edges such 
a clump will trap photons for many sweeps as their travelled 
distance after N sweeps is given by y/ND where T) is small. 
This is another reason against the use of high sampling powers 
throughout the simulation domain. 



7. Discussion 

Now that we have identified and corrected the various effects due 
to anisotropy in the Delaunay grid we can take a step back and 
discuss their implications in a broader context. 

7. 1 . Prevention versus correction 

It is impossible to define the perfect grid in general. As already 
pointed out, a large dynamic range inevitably leads to gradients 
that, in turn, give rise to the unphysical effects described in this 
text. Although we showed that these effects can be corrected for 
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by weighting schemes, these corrections are not fail-safe if the 
gradients are too strong. In particular, when the density of nu- 
clei changes appreciatively on length scales smaller than T> (or, 
more precise, Q„ < 1), the lack of edges towards the under- 
dense region may lead to extreme decollimation as there simply 
are hardly any edges that point that way. In this case, a weighting 
scheme is not feasible as there are no edges to assign weights to 
in the under-dense direction. Simply said, we should start with a 
fairly well behaved grid in the first place. 

One could argue that if the grid is constructed such that 
the unphysical effects are below a predefined tolerance, the de- 
ployment of weighting schemes can be circumvented. On the 
other hand, this puts rather stringent constraints on the grid and 
it might be preferable to retain more dynamic range (and thus 
structure) in the grid at the expense of the need for weighting. 

An intermediated solution will often be the most viable op- 
tion. Retaining the dynamic range of the data as much as possi- 
ble while weighting schemes and DCT are employed at certain 
locations in the simulation domain only. 

7.2. Number of steps versus L eq and L def 

In Sec. [3] we have derived typical length scales at which the un- 
physical effects become important in a relative sense. Being cau- 
tious, we have set these length scales to many times the size of 
the simulation domain but this is not necessary in most real life 
situations. This is because radiation typically travels only moder- 
ate distances along the grid. It might therefore be more intuitive 
to think about the number of steps a photon travels along the grid 
before being scattered or absorbed instead of L eq and Ldef ■ In ei- 
ther the event of absorption or scattering, the past of the photon 
is essentially erased and it starts a new life, so decollimation, 
deflection and drift are basically 'reset' to zero. 

We can thus assess the possible harm of anisotropy in a sim- 
ulation by monitoring the fraction of Delaunay length and phys- 
ical mean free path, R, defined as 



R = 



D 



(41) 



mfp 



First of all, R should not be much smaller than 1, because 
this would result in decollimation and in extreme cases even in 
simulating straight light rays by a random walk (see also the 
discussion in Sec. |3.3.lj ). 



In the limit of large R, however, a large portion of the photons 
will be absorbed after travelling along this edge. This is merely 
a resolution issue; the photons will never travel further than one 
Voronoi region away from the place where they ought to have 
been absorbed. 



Using a sampling function of the form of Eq.(37 i would re- 
sult in R being a constant along the grid. If we use a more gen- 
eral sampling method, R is no longer a global constant. We can 
however still adjust the local optical depth for absorption with 
the local value of R. If we do so, the place where a photon fi- 
nally gets absorbed or leaves the simulation domain still has the 
right probability distribution, but the interpretation of interme- 
diate photon densities in the simulation becomes less intuitive. 
Already the photon densities are no 'time' densities but 'step' 
densities, but now photons in regions with different densities will 
travel a different amount of mean free paths per sweep. Still, if 
we are interested in the effect of the photons on the medium and 
in the directions of photons that leave the box, we are justified 
to sample with any other power than 3, as long as we make local 
corrections to the optical depth proportional to R. 



7.3. A dynamic grid 

We have argued that the ratio between the local mean free path 
and the Delaunay edges should be chosen with care in order to 
ensure meaningful results. Because the local mean free path can 
change due to changes in the medium properties, the Delaunay 
edge length should adapt accordingly to keep R fixed. 

Recently the SimpleX code has been extended to handle dy- 
namic media (Paardekoo per et al.|2009) , meaning that nuclei can 
be added or deleted during the simulation. This allows the trian- 
gulation to adjust itself to the (changing) local properties of the 
medium. 

We can set up strict conditions under which a nucleus gets 
deleted (and its attributes divided between its neighbours). A 
possible condition could be: delete a nucleus whenever the lo- 
cal mean free path of the ionising photons exceeds the average 
length of its outgoing edges. In this way we keep deflection and 
decollimation acceptably small in the case of ballistic transport. 
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8. Future work 

In this section we discuss several limitations of the current 
SimpleX method and indicate a possible route to overcome these 
limitation in future versions of the code. 



8.1. Multiple frequencies 

In the current incarnation of SimpleX, a grey approximation is 
used where the effect of a source spectrum is taken into ac- 
count by the use of effective, intensity weighted, cross sections. 
Implementation of a full multi-frequency treatment is not triv- 
ial as it would formally require a separate grid for every new 
frequency bin. In practice, however, we have found that by 
transporting similar frequencies over the same grid, the com- 
putational demands will not become prohibitive. A version of 
SimpleX incorporating multiple frequencies is currently under 
development and we postpone a discussion of related issues to 
future articles. 



8.2. Cosmological redshift 

For cosmological simulations, the effect of expanding space will 
result in shifting of frequency of the transported photons. In 
our current one-frequency approach this is nothing more than a 
gradual loss of ionizing energy for the photons. In future multi- 
frequency versions of the code more care should be taken of this 
issue especially when line-transfer is included. 

In principle, also relativistic 'beaming' could play a non- 
neglegible role in simulations of large cosmological volumes. 
Due to expansion of the universe, an absorbing parcel of gas 
will have different relative velocities with sources at different 
distances. The resulting enhancement of the observed intensity 
(e.g. Eq. 4.9 in Rybicki & Lightman 1986) can be substantial. 
Such beaming could in principle be modeled in SimpleX by ad- 
justing the weighting factors for diffuse and ballistic transport 
appropriately but we leave a detailed analysis for future work. 



tematic effects. Although only some of the described effects 
(diffuse drift and ballistic decollimation) need to be consid- 
ered in typical applications (see Paardekoo per et aL| [2009 



8.3. Strong scattering 

Finally, in the case of strong scattering (if the cross section for 
absorption is negligibly small compared to the cross section for 
scattering) the SimpleX method only converges asymptotically to 
the correct answer. Colloquially speaking this is due to the fact 
that transport happens on a cell-to-cell basis and this can take 
a large number of sweeps to bring all scattering particles in the 
simulation domain into equilibrium. This problem is well known 
and has been overcome by classical radiative transfer methods 
by the use of so-called accelerated Lambda iteration methods. 

In the cosmological applications of SimpleX, strong scatter- 
ing is not a dominant process, however, and we will not pursue 
an implementation of accelerated Lambda iteration schemes in 
our method in the near future. 



9. Summary 

1 . The SimpleX algorithm is a relatively new and non-standard 
method and as such deserves a thorough analysis of its sys- 
tematic effects. This assessment has two important purposes: 
First, exploring the boundaries of its reliability and accuracy 
and, second, increasing its region of application through im- 
provements. 

2. The mathematically transparent nature of the SimpleX algo- 
rithm allows for a rigorous and general assessment of its sys- 



we have investigated all four in detail for completeness. 

3. The use of a random Delaunay triangulation in the radia- 
tive transfer method SimpleX introduces global errors when 
the point-density is inhomogeneous. We have identified and 
quantified four distinct effects: 'drift' and 'clustering' of 
photons in diffusive transport, and 'decollimation' and 'de- 
flection' in ballistic transport. 

4. We show that decollimation become troublesome only when 
the number of traversed edges becomes larger than roughly 
5 steps. This implies that one either preclude this regime or, 
when this is undesirable or impossible, correct for the un- 
physical behaviour. 

5. A weighting method for ballistic transport has been shown 
to decrease the effect of decollimation significantly pushing 
the limit for correct transport from 5 to 12 steps. 

6. We have shown how drift and clustering can be adequately 
corrected for by adopting a weighting scheme. Three weight- 
ing schemes for diffusive transport have been discussed and 
compared. 

7. The use of a direction conserving variant of the ballistic 
transport has been introduced as an alternative to the weight- 
ing procedure and to prevent loss of angular direction in op- 
tically thin regions. This approach provides a means of cor- 
recting for deflection, drift and clustering as well. The com- 
putational cost and need for additional randomisation of the 
solid angle directions makes the employment of DCT not 
generally favourable, however. 

8. Finally, we have applied our analysis to an astrophysically 
relevant example. This elucidates the intimate relation be- 
tween the construction of a computational grid and the pos- 
sible need for weighting schemes. 
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